In [1]:
# Setting up a model and a mesh for the MT forward problem

In [24]:
import SimPEG as simpeg
from SimPEG import MT
sys.path.append('/home/gudni/gitCodes/python/telluricpy')
import telluricpy

In [6]:
# Define the area of interest
bw, be = 556500, 558000
bs, bn = 7133500, 7133500
bb, bt = 0,480

In [14]:
# Build the mesh
# Design the tensors
hSize,vSize =  50., 12.5
nrCcore = [10, 8, 6, 4, 2, 2, 2, 2, 2]
hPad = simpeg.Utils.meshTensor([(hSize,9,1.5)])
hx = np.concatenate((hPad[::-1],np.ones(((be-bw)/hSize,))*hSize,hPad))
hy = np.concatenate((hPad[::-1],np.ones(((bn-bs)/hSize,))*hSize,hPad))
airPad = simpeg.Utils.meshTensor([(vSize,13,1.5)])
vCore = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcore,(simpeg.Utils.meshTensor([(vSize,1),(vSize,8,1.3)])))])[::-1]
botPad = simpeg.Utils.meshTensor([(vCore[0],8,-1.5)])
hz = np.concatenate((botPad,vCore,airPad))
# Calculate the x0 point
x0 = np.array([bw-np.sum(hPad),bs-np.sum(hPad),bt-np.sum(vCore)-np.sum(botPad)])
# Make the mesh
meshFor = simpeg.Mesh.TensorMesh([hx,hy,hz],x0)

In [17]:
print np.sum(vCore)
print meshFor.nC
print meshFor


1137.29994775
164256
  ---- 3-D TensorMesh ----  
   x0: 550883.50
   y0: 7127383.50
   z0: -8191.26
  nCx: 48
  nCy: 58
  nCz: 59
   hx: 1922.17, 1281.45, 854.30, 569.53, 379.69, 253.12, 168.75, 112.50, 75.00, 30*50.00, 75.00, 112.50, 168.75, 253.12, 379.69, 569.53, 854.30, 1281.45, 1922.17
   hy: 1922.17, 1281.45, 854.30, 569.53, 379.69, 253.12, 168.75, 112.50, 75.00, 40*50.00, 75.00, 112.50, 168.75, 253.12, 379.69, 569.53, 854.30, 1281.45, 1922.17
   hz: 2613.29, 1742.19, 1161.46, 774.31, 516.20, 344.14, 229.42, 152.95, 2*101.97, 2*78.44, 2*60.34, 2*46.41, 2*35.70, 4*27.46, 6*21.13, 8*16.25, 10*12.50, 18.75, 28.12, 42.19, 63.28, 94.92, 142.38, 213.57, 320.36, 480.54, 720.81, 1081.22, 1621.83, 2432.74

In [35]:
# Save the mesh
meshFor.writeVTK('nsmesh.vtr',{'id':np.arange(meshFor.nC)})
nsvtr = telluricpy.vtkTools.io.readVTRFile('nsmesh.vtr')

In [ ]:
topoSurf = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile('../Geological_model/CDED_Lake_Coarse.vtp'))
activeMod = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(nsvtr,topoSurf)

In [39]:
telluricpy.vtkTools.io.writeVTUFile('activeModel.vtu',activeMod)

In [40]:
# Get active indieces 
activeInd = telluricpy.vtkTools.dataset.getDataArray(activeMod,'id')

In [83]:
# Make the conductivity dictionary
# Note: using the background value for the till, since the extraction gets the ind's below the till surface
geoStructFileDict = {'Till':1e-4,
 'XVK':3e-2,
 'PK1':5e-2,
 'PK2':1e-2,
 'PK3':1e-2,
 'HK1':1e-3,
  'VK':5e-3}

In [75]:
# Loop through
extP = '../Geological_model/'
geoStructIndDict = {}
for key, val in geoStructFileDict.iteritems():
    geoPoly = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile(extP+key+'.vtp'))
    modStruct = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(activeMod,geoPoly,extBoundaryCells=True,extInside=True,extractBounds=True)
    geoStructIndDict[key] = telluricpy.vtkTools.dataset.getDataArray(modStruct,'id')

In [89]:
# Make the physical prop
sigma = np.ones(meshFor.nC)*1e-8
sigma[activeInd] = 1e-3 # 1e-4 is the background and 1e-3 is the till value
# Add the structure
for key in ['Till','XVK','PK1','PK2','PK3','HK1','VK']:
    sigma[geoStructIndDict[key]] = geoStructFileDict[key]

In [90]:
# Save the model
meshFor.writeVTK('nsmesh_0.vtr',{'S/m':sigma})

In [91]:
# Set up the forward modeling
freq = np.logspace(5,0,26)
np.save('MTfrequencies',freq)

In [58]:
# Find the locations on the surface of the model.
# Get the outer shell of the model
actModVTP = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.extraction.geometryFilt(activeMod))
polyBox = vtk.vtkCubeSource()
polyBox.SetBounds(bw,be,bs,bn,bb,bt)
polyBox.Update()
# Exract the topo of the model
modTopoVTP = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(actModVTP,telluricpy.vtkTools.polydata.normFilter(polyBox.GetOutput()),extractBounds=True)

In [92]:
# Make the rxLocations file
x,y = np.meshgrid(np.arange(bw+25.,be,50),np.arange(bs+25.,bn,50))
xy = np.hstack((x.reshape(-1,1),y.reshape(-1,1)))
# Find the location array
locArr = telluricpy.modelTools.surfaceIntersect.findZofXYOnPolydata(xy,modTopoVTP)
np.save('MTlocations',locArr)

In [60]:
# Running the forward modelling on the Cluster.
# Define the forward run in findDiam_MTforward.py

In [95]:
%matplotlib qt
sys.path.append('/home/gudni/Dropbox/code/python/MTview/')
import interactivePlotFunctions as iPf

In [99]:
# Load the data
mtData = np.load('MTdataStArr_nsmesh_0.npy')
mtData


Out[99]:
array([ (100000.0, 556525.0, 7133025.0, 442.5, (-0.4218512292686255+0.3823156030050135j), (-46.92708407644145-28.827243176833832j), (49.954240237816684+28.304900026316226j), (0.6587762718014313-0.3911370761200063j), (-0.003324532589959883-0.008319925856157012j), (-0.01850717106625786+0.020712450433413677j)),
       (100000.0, 556575.0, 7133025.0, 442.5, (-0.01726045031745334+0.1943430384592769j), (-46.74238696054156-28.742861178661375j), (50.071153918121794+28.361729613106952j), (-0.03927194461508535-0.24754114198717278j), (-0.007264019653180863-0.007530857913013391j), (-0.02104071911847556+0.022778095807779033j)),
       (100000.0, 556625.0, 7133025.0, 442.5, (0.3715278792570291+0.007127420821101053j), (-46.94950392918378-28.6043071241115j), (49.82051013844956+28.62941149137527j), (-0.6718553350812287-0.11045535419742107j), (-0.010420762937427592-0.008463498923996504j), (-0.01871299493368591+0.02405049231561648j)),
       ...,
       (1.0, 557875.0, 7134975.0, 442.5, (0.5425667495297678+0.19935812310850504j), (-2.026685596578978-0.7604536835983275j), (1.7111876926366725+0.6357136615294343j), (-0.5962284716699029-0.22220234677435738j), (0.0020321487290150837+0.0007776291784596612j), (0.0006279739039880312+0.00026577865648246327j)),
       (1.0, 557925.0, 7134975.0, 442.5, (0.39234841712090124+0.14459524436646765j), (-1.5259801336396221-0.5747168781273916j), (1.5103939959937822+0.5625349703912992j), (-0.44931421262228693-0.16783373592862247j), (0.0041757731387068045+0.0015798277141383437j), (0.002809153001429798+0.001090843492842591j)),
       (1.0, 557975.0, 7134975.0, 455.0, (-0.15265548345520036-0.05718918005767595j), (-1.0044265135244046-0.3780513949684229j), (1.043479913030608+0.3879577740540316j), (0.138555232654058+0.05252180227205386j), (0.004572389599954534+0.0017197152959960568j), (0.00304885579854783+0.0011763896890463455j))], 
      dtype=[('freq', '<f8'), ('x', '<f8'), ('y', '<f8'), ('z', '<f8'), ('zxx', '<c16'), ('zxy', '<c16'), ('zyx', '<c16'), ('zyy', '<c16'), ('tzx', '<c16'), ('tzy', '<c16')])

In [110]:
iPf.MTinteractiveMap([mtData])


Out[110]:
<interactivePlotFunctions.MTinteractiveMap at 0x7f49d3bc16d0>

In [104]:
# Looking at the data shows that data below 100Hz is affected by the boundary conditions, 
# which makes sense for very conductive conditions as we have.
# Invert data in the 1e5-1e2 range.

Run the inversion on the cluster using the inv3d/run1/findDiam_inversion.py


In [106]:
drecAll = np.load('MTdataStArr_nsmesh_0.npy')

In [109]:
np.unique(drecAll['freq'])[10::]


Out[109]:
array([    100.        ,     158.48931925,     251.18864315,
           398.10717055,     630.95734448,    1000.        ,
          1584.89319246,    2511.88643151,    3981.07170553,
          6309.5734448 ,   10000.        ,   15848.93192461,
         25118.8643151 ,   39810.71705535,   63095.73444802,  100000.        ])

In [ ]: